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ABSTRACT 


Starting with a least squares formulation of the parameter estimation problem, 
both fixed data and data-adaptive iterative algorithms are developed. We apply two 
new techniques, namely diagonal pertubation and multiple partitioning, to existing 
finite impulse response (FIR) and infinite impulse response (IIR) fixed data matrix 
splitting algorithms, resulting in improved performance. Also, we extend the fixed 
data algorithms to the data-adaptive case, and contrast them with FIR and IIR 
recursive least squares (RLS) algorithms. Computer simulations are used to evaluate 
the computational effectiveness of the new algorithms. We show the general rate 
of convergence for the algorithms, evaluate their ability to correctly represent the 
spectral components of simulated system frequency response in noise, and present 
system performance, when the order of the model is chosen to be larger than the 


known system order (over-modeling). 
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I. INTRODUCTION 


The general field of system modeling and parameter estimation is a rich and 
current area of research. In this thesis we apply four principal concepts to solving the 


least squares normal equations, which have the general form: 
Ra=r (1.1) 


where RF is the data correlation matrix, a is the parameter vector to be determined, 
and r is the cross-correlation vector. First, we desire to iteratively solve the problem 
rather than use direct inversion of the correlation matrix R. Second, we use a Toeplitz 
approximation matrix splitting technique [Ref. 1, 2] to set up the iterative equations. 
Third, we increase the diagonal dominance of the correlation matrix R in order to 
improve the convergence properties of the iterative equations. Finally, we partition 
the normal equations (1.1) in order to improve the computational efficiency and rate 
of convergence of the iterative algorithms. In all these cases our interest is to study 


both fixed data or off-line algorithms and data-adaptive or online algorithms. 


A. PERFORMANCE CRITERIA 

The performance criterion for the algorithms developed in this thesis is based 
on a least squares formulation. The algorithms minimize the sum of the squared 
error between the true system output and the estimated system output. We choose 
to use the covariance method for setting up the problem in order to remain within the 
available data. Although the covariance method causes the data correlation matrix 


R of (1.1) to be non-Toeplitz, it has the potential to provide a less biased parameter 


estimate than the autocorrelation method [Ref. 3]. 





B. THESIS OVERVIEW 

This thesis is divided into five chapters, including the Introduction. In Chapter 
II we explicitly develop the least squares data formulation for both the finite impulse 
response (FIR) and infinite impulse response (IIR) system models. For the IIR system 
we take special care to account for the dependency of the filter output on the input 
and output filter coefficients. Chapter III presents iterative algorithms for FIR sys- 
tems. First we consider the fixed data case, and present the simplistic Gauss—Siedel 
algorithm as a way of introducing the concepts of matrix splitting and iteration. Then 
we present the Toeplitz approximation matrix splittiug algorithm, and develop two 
modified versions of it. Here we apply two new techniques: diagonal pertubation and 
multiple partitioning. Next, we develop and contrast a new data-adaptive algorithm 
with the well known recursive least squares (RLS) algorithm. IIR systems are inves- 
tigated in Chapter IV. First we present the fixed data IIR Toeplitz approximation 
iterative algorithm, and then attempt to apply diagonal pertubation and partition- 
ing to it. Achieving marginal success for the fixed data case, we then turn to the 
data-adaptive case. Here we develop a new data-adaptive IIR algorithm based on the 
Toeplitz approximation matrix splitting that was used in the fixed data case. This 
algorithm performs well, and does not have the stability problems associated with 
the IIR RLS algorithm, which we present for comparison. In the final chapter, we 


summarize the results of simulation and recommend topics for future research. 


II. LEAST SQUARES DATA FORMULATION 


The central problem of this thesis is to determine the parameters of an optimal 
linear filter based on the input and output measurements of a system. For the system 
shown in Figure 2.1, both the input data sequence z(n) and the output data sequence 
y(n) are known. The objective is to determine a linear filter, or model system which 
produces the output from the given input. Using least squares minimization tech- 
niques gives a filter that is optimal when it minimizes the sum of the squared error 
between the known output y(n) and the filter output y(n). The first section of this 
chapter contains the least squares data formulation for the finite impulse response 
(FIR), or non-recursive filter. The next section presents the least squares data for- 
mulation for the infinite impulse response (IIR), or recursive filter, which has some 
subtle but significant differences from the FIR case. This chapter concludes with 


some generalizations about the algorithms to be developed. 







MODEL 





SYSTEM 


x(n) ERROR 


UNKNOWN 






y(n) 





SYSTEM 


Figure 2.1: Fundamental System Model 


A. FINITE IMPULSE RESPONSE (FIR) MODELS 
Consider an M**-order FIR filter with input z(n) and output y(n). The output 


at time n is simply a weighted linear combination of the input: 
j(n) = x(n)oM 4 x(n —1)a + ---+2(n— Mah. (2.1) 


This can be written as 
§(n) = xna™, (2.2) 


where x, = (x(n) z(n—1) --- 2(n~M)] isa 1x M+1 input data vector, aM@ = [aM aM 


--- aM]? is a M+1x 1 vector containing the filter weights, or coefficients, and 9(n) is 
the filter output at time n. The superscript M indicates that these are the coefficients 
of an Mt*-order filter. Because the true value of a” is not known, equation (2.2) 
produces an estimate for the output, y(n). Consequently, the error e(n) between the 


true output y(n) and the FIR filter output y(n), is given by 
e(n) = y(n) — §(n) = y(n) — xna™. (2.3) 


The key to a least squares solution of equation (2.2) for a“ is forming an 


overdetermined set of P + 1 equations (where P > M is necessary for a unique 


solution): 
y(n) z(n) z(n—1) ate z(n ~ M) aM 
= g(n — 1) ia z(n-1) 2(n-2)  --- 2(n—-M-—1) aM 
j(n — P) en =P) w= P =) eo z(n — M — P) all 
Xn 
ae ete a (2.4) 
Xn-P 


this can be written more compactly in matrix notation as 





Equation (2.5) is overdetermined in the sense that there are more equations than 
there are unknowns, or, equivalently, there are more rows of the data matrix X than 
there are coefficients in the filter coefficient vector aM. Using (2.5), the error vector 
can be written as: e = y—y = y—Xa™. An optimal least squares solution minimizes 


the sum of the squared errors ¢, which is represented as: 


n 


C= DO lei? = e7e 
is (y — Xa”) ly — Xa") 
(y7 — a? XT)(y — Xa”) 
yTy —y™Xa —aM@7xTy4aM7xTXa. (2.6) 


The minimization is accomplished by taking the derivative of the sum of the squared 
errors ¢ with respect to the filter coefficient vector a” [Ref. 4] and setting it equal 
to zero: 


Bet) = 0—2XTy + 2X7Xa™ =0. (2.7) 


Rearranging gives 

XTXaM = XTy (2.8) 
as the requisite condition for ¢ to be minimum. It is important to notice that the input 
data which forms X, and the true output data contained in y are independent of the 
filter coefficients contained in a“. Because (2.8) contains the estimate of the optimal 
M*-order FIR filter coefficients a”, solving (2.8) determines the optimal FIR filter. 
Therefore, equation (2.8) defines the FIR filter problem to be solved in Chapter III. 


A more standard representation of (2.8), referred to as the normal equation is 
Ra” =r, (2.9) 


where R = X7X is referred to as the correlation matrix, and r = X‘y is the cross- 


correlation of the input data matrix X and the output data vector y. Theoretically 


5 


if R is full rank, the solution to (2.9) is given by a” = R-!r. This thesis presents 
iterative algorithms which are alternatives to the direct inversion of R. Directly 


inverting R to solve (2.9) has three major shortcomings: 


1. Direct inversion is computationally intensive; on the order of M* multiplications 


for a M x M correlation matrix R. 


2. If R is nearly singular, or very poorly conditioned it may not be possible to 


compute Ro}. 


3. Unlike iterative algorithms, direct inversion cannot provide an incremental es- 


timate of the solution. 


B. INFINITE IMPULSE RESPONSE (IIR) MODELS 

The general causal linear IIR filter is more complicated than the FIR filter 
presented above. This is because the filter output y(n) is a linear combination of the 
input z(n), z(n —1),---, z(n — M), as well as the previous output 9(n — 1), y(n — 2), 
-++, 9(n — N), where M and N are the order of the input and output coefficients, 


respectively. Starting with the difference equation representation of the IIR filter: 


Gn) = Gn —1)by + G(r — 2)by + --- +. 9(n — N)bH 
+ 2(n)ai’ + 2(n—1)a! +---4+2(n- M)a™, (2.10) 


b” 
in) = [gna : wl] 5 


y(n) = 28, (2.11) 


which can be written as 


where Z, = [Yn-1 : Xn] = [9(n — 1) #(n —2)--- g(n—N) : 2(n) z(n—1)---2(n—M)] 
[ 


isa x N+M +1 data vector, 9 =(bN7 : aM7]T = [BN WN...oN i gM gM... MJT 


6 





isa N+ M+1 x 1 vector containing the IIR filter weights. Again, because the true 
value of b™ and a™ are not known, equation (2.11) will only produce an estimate 
for the output, j(n). Thus, the error between the true output y(n) and the IIR filter 


output g(r), is given by 
e(n) = y(n) —J(n) = y(n) — 29. (2.12) 


As before, the least squares solution to (2.11) for @ comes from an overdeter- 


mined set of P + 1 equations (where P > N + M is necessary for a unique solution): 


y(n) 
F y(n — 1) 
y= ; 
y(n — P) 
by 
Hm=1) MM) atm) atm) | 
_ y(n—-2) +++ #(mn—-N-1) : a2(n-1) --- 2(n—M-—1) xy 
: : : : aM 
j(n-—P-1) --- 9j(n-N-P) : z(n-—P) --- e(n—M-—P) : 
aM 
M 
Yn-1 : Xn Zn 
7 : Zn- 
PS OTR: Fg a ley (2.13) 
Yn-P-1 : Xn—P eae 
or, in matrix notation 
y= 20, | (2.14) 


where the data matrix Z = [27 27_, --- 27_p]?, which is used to form an error vector, 
e = y~ Z@. Once again, the optimal least squares solution minimizes the sum of the 


squared errors ¢, which now has the form: 


C= Yo eG)? = ee 

(y — 26)" (y — Ze) 

(y7 — @72Z7)(y ~ Ze) 

yTy —y™Zo —oTZTy +0727 Ze. (2.15) 


lI 


If the gradient terms that arise from formal differentiation are neglected, the derivative 


of the sum of the squared errors ¢ with respect to the filter coefficient vector @ is: 


WO) _ ) _ o9T Top 
56) 7° QZTy +2Z7 Ze =0. (2.16) 


Rearranging gives 


Z'Zea=Z'y (2.17) 


as the requisite condition for ¢ to be minimum. Again, (2.17) contains the estimate 
of the IIR filter coefficients 6, and its solution determines the optimal least squares 
IIR filter. Therefore, equation (2.17) defines the IIR filter problem to be solved in 


Chapter IV, and the normal equation representation is 
Re =r, (2.18) 


where R = Z?Z is referred to as the correlation matrix, and r = Z’y is the cross- 
correlation vector. 

The differentiation in (2.16) appears quite simple because it ignores the de- 
pendence of the filter output y(n) on the previous filter input and output. Formal 


differentiation of ¢ gives [Ref. 5]: 


BY Djon-P(y(3) — 9(5))? 
$i «| Perc Tide, Sian 
(| bu |) Lal E3-.-r(vG) 36)? 





Letting 6, = aye) and a, = a) represent the vector partial derivatives of tne filter 


output, we can write 
n 


De | i | (y(3) — 236) = 0, (2.20) 


j=n— 


or, rearranging and letting #; = [67 a7]” gives 


» wm] = > by(3) 
j=n—P j=n-P 


J 


a = vy, (2.21) 
r 


so that we have the same form as the normal equations of (2.18), Re =r. However, 
R no longer represents the autocorrelation of the data matrix Z. In this case R 
represents the cross-correlation matrix between the gradient components of the filter 
output W and the filter data matrix Z, and r represents the cross-correlation between 
WY and the true output vector y. 

Regardless of the formulation, a standard solution to (2.18) is given by @ = 
R-'r. Again, the objective is to present iterative algorithms which are alternatives 
to the direct inversion of R. Finally, for the algorithms presented in this thesis, there 


are three general characteristics governing their suitability: 


1. The algorithm should use the minimum amount of data required to produce a 


solution. 


2. The algorithm should converge quickly. The desired number of iterations of the 
algorithm should not exceed 2 to 3 times the number of parameters being solved 


for. 


3. The computational complexity of the algorithm should be less than that of 


direct inversion. 








Ill. FINITE IMPULSE RESPONSE (FIR) 
SYSTEMS 


In this chapter we consider the finite impulse response (FIR) system. We develop 
two general categories of solution algorithms, fixed data and data-adaptive. The fixed 
data algorithms are the Gauss-Seidel iterative method and the Toeplitz approxima- 
tion iterative method. The data-adaptive algorithms are the recursive least squares 
(RLS) method and Toeplitz approximation method. Before developing the solution 
algorithms, it is worth noting that the basic problem formulation described in Chap- 
ter II is a covariance rather than a correlation formulation of the data; therefore, the 
matrix Rin Ra! = r is symmetric but not Toeplitz. Accordingly, we cannot apply 
the Levinson recursion algorithm (Ref. 6] which requires Toeplitz structure in the 
correlation matrix. In our formulation, the data matrix is formed using the covari- 
ance method, which has the potential to provide a less biased least squares solution 


than a correlation method formulation [Ref. 3]. 


A. FIXED DATA ALGORITHMS 
1. Gauss-Seidel Method 
A very simple and straightforward iterative algorithm is the Gauss~Seidel 
method [Ref. 7]. We drop the superscript M from a™ for simplicity. Unless otherwise 
stated, we are considering an M‘*-order FIR system. Starting with (2.9), Ra =r, 
split R into L ++ D+U, where L is a matrix containing the strictly lower triangular 
elements of R, U is a matrix containing the strictly upper triangular elements of R, 


and D is a diagonal matrix containing the main diagonal elements of R. Substituting 
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this into (2.9) gives 
(L+D+U)a=r, (3.1) 


or 


(L+D)a=—Uatr. (3.2) 


Making (3.2) into an iterative algorithm requires that we isolate the parameter vector 
a on the left side of the equation. Pre-multiplying both sides of the equation by 
(L+D)7}: 

(L+ D)"(L+ D)a=(L+D)"'(-Ua+r) (3.3) 


we have 
a=—(L+D)'Ua+(L+D)"'r. (3.4) 
Note that (Z + D) must be nonsinglar for this technique to work. Now there are two 
occurrences of the parameter vector a. We define the one on the right side of (3.4) as 
the current value a‘), and the one on the left side as the updated value a‘*+), where 
k is the iteration index. Thus, the final form of the Gauss-Seidel iterative algorithm 
is: 
al*+) = —(L + D)-1Ual® + (L+ D)"r. (3.5) 
Unfortunately, this algorithm takes an excessive number of iterations to converge to 
the true parameter values when F is even moderately ill-conditioned. In an effort 
to improve this performance, we applied the successive overrelaxation (SOR) accel- 
eration technique (Ref. 9] to (3.5). The SOR technique requires that we add an 


acceleration parameter w as follows: 
al®*)) = —(wL + D){(1 —w)D — wa + (wh + D)“ wr. (3.6) 


Notice that when w equals 1, (3.6) degenerates into the standard Gauss-Seidel algo- 


rithm. In order to speed up the rate of convergence of the algorithm, the acceleration 
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parameter w must assume its optimum value, which is given by [Ref. 9]: 


2(1 = yi = priae) = 
Wopt = i wa pS (3. 7) 
Umar 
In (3.7) the term maz represents the largest eigenvalue of the Jacobi matrix, which 
has the form D~!(—ZL ~— U). Unfortunately, though this is a popular and usua!ly 
successful acceleration technique, for the basic problem of this thesis it was found to 
produce little improvement in the performance of the Gauss-Seidel algorithm. 


In order to observe the performance of the algorithms developed in this 


chapter, they are simulated with the following FIR system: 


y(n) = 0.52(n —1) + 0.252(n — 2) ~ 0.52(n — 4) + 0.0532(n — 5) + 0.03452(n — 6) 
— 0.76x(n — 7) + 3.52(n — 8) ~ 1.0032z(n ~ 9) — 0.00312(n — 10). (3.8) 


Figure 3.1 shows the generally slow and gradual rate of convergence of the Gauss- 
Siedel and Gauss-Siedel with SOR algorithms. Each parameter track line of Figure 
3.1 corresponds to a coefficient in (3.8). Notice that as the number of iterations k 
increases, the parameter track lines flatten out and approach the values in (3.8). 
2. Toeplitz Approximation Algorithm 

One reason for the slow rate of convergence of the Gauss-Seidel method is 
that it fails to capitalize on the structure of the covariance matrix R. The covariance 
matrix is symmetric but not Toeplitz. However, it does have an underlying Toeplitz 
structure. This point becomes clear when the number of data points used to form 
the matrix is increased; the covariance matrix then approaches a Toeplitz structure. 
Toeplitz matricies have many nice properties, and efficient methods such as the Levin- 
son recursion are available to invert a Toeplitz matrix. An iterative algorithm which 
takes advantage of the near Toeplitz structure of the covariance matrix [Ref. 1] is de- 


veloped by splitting R into a Toeplitz matrix T, and a residual matrix S. The matrix 
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Figure 3.1: Parameter tracks and general rate of convergence for the stan- 
dard FIR Gauss-Seidel and Gauss-Seidel with SOR algorithms. 


T is obtained by averaging the diagonal elements of R. This Toeplitz approximation 
provides a natural splitting of the covariance matrix, R = T + S, which is used to 


develop an iterative algorithm as follows. Beginning with 


Ra=r, (3.9) 
substituting for R = T + S gives 
(T+S)a=r, (3.10) 
or, rearranging we have, 
Ta=r-—Sa. (3.11) 
But S = R—T, thus 
Ta=r—-(R-T)a (3.12) 
and isolating a gives 
a=T'r-T'Ra+a. (3.13) 
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Finally, the iterative algorithm las the form: 
al*t}) = Typ — T71Ral®) 4 al), (3.14) 


And if we let a9 = T~'r be the initial estimate of the parameter vector a, then we 


have: 


a+) Say+al)—T- Ral), (3.15) 


The Toeplitz approximation matrix T can also be viewed in another light 
(Ref. 8]. It can be thought of as a pre-conditioning matrix for R. Indeed, part 
of the difficulty with the Gauss-Seidel method of the previous section is that as R 
becomes more ill-conditioned, the iterative algorithm takes 1. uch longer to converge. 
We observed that the product T-!R does in fact have a bette: condition number than 
R by at least an order of magnitude when R is very ill-conditioned. Thus, we might 
expect that the algorithm of (3.15) will converge to the true parameter vector a faster 
than (3.5) or (3.6). This is in fact the case. In Figure 3.2 it is clear that the Toeplitz 
approximation algorithm gives much faster convergence to the true parameter values 
than the Gauss-Seidel algorithm. 

3. Toeplitz Appreximation with Diagonal Pertubation Algorithm 

One significant shortcoming of (3.15) is that it fails to converge as the 
covariance matrix R becomes extremely ill-conditioned. Despite the fact that T 
functions as a pre-conditioner for R, the algorithm of (3.15) will not converge as the 
spectral radius (or largest eigenvalue) of (T~'S) becomes greater than one [Ref. 9]. 
Generally, the condition number of & is an indication of 1ow poorly conditioned it 
is. We observe that as the number of data points used to form R is reduced from 
a large number (i.e., 5 to 10 times the order of the filter being solved for) to near 
the minimum number required for a unique solution (i.e., the order of the filter), the 


condition number of # grows quite large. Under these circumstances the algorithm of 


14 


Gauss-Siedel, P=50 





Parameter Value 
Parameter Value 





0 5 10 15 20 


Iterations, k Iterations, k 


Figure 3.2: Parameter tracks and general rate of convergence for the FIR 
Toeplitz approximation and the Gauss-Seidel algorithms. 

(3.15) does not converge. By constrast the Gauss-Seidel algorithm always converges 
to the true parameter values, although it takes an excessivly large number of iterations 
for the parameter values to approach their true values (see Figure 3.3). Viewing the 
ill-conditioned nature of R as an obstacle to satisfactory performance of (3.15), we 
now introduce a technique which improves the condition of R and permits a modified 
version of (3.15) to converge as the number of data points approaches the minimum 
required. 

Beginning with the fact that any identity matrix J is invertible with a 
condition number of one, we seek to alter R such that it becomes more like an identity 
matrix. This is accomplished by adding a diagonal matrix to R such that the resulting 
matrix has its main diagonal elements at least an order of magnitude larger than any 
other element in the matrix, then this new matrix has increased diagonal dominance 
and begins to look more like an identity matrix. The consequence of applying this 


technique is that we improve the condition of R and the performance of (3.15). 


15 


Gauss-Siedel, P=14 





3 
g g 
a a 
> > 
3 5 
a 3 
A. a 
Iterations, k Iterations, k 


Figure 3.3: Parameter tracks and general rate of convergence for the FIR 
Toeplitz approximation and Gauss-Seidel algorithms. 


If we let o be a scalar multiplier and the Toeplitz approximation splitting 
of & be defined as R = 7, +S, then we can derive a modified version of (3.15) as 
follows. Beginning with 


Ra=r, (3.16) 


adding a diagonal matrix (scaled by o) to R gives 
(R+olha=r+oa (3.17) 
and substituting for R = T, + S we have 
(T.+S+oelja=r+ea. (3.18) 
Now, letting T = T, + oJ, we can write 
Ta=r+oa—Sa (3.19) 
but S= R+ol—T, and substituting gives 


Ta=r+oa-—(R+ol—T)a. (3.20) 
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After rearranging it is easy to see that the oa terms cancel 
Ta=r+oca—oa—(R-Tha (3.21) 


so that we can now write 


a=T'r—T'Ra+a. (3.22) 


Now, letting ag = T—'r and rearranging, the final form of the Toeplitz approximation 


with diagonal pertubation iterative algorithm is: 
al*tl) = ay + al®) — T-Ral®), (3.23) 


Interestingly, the algorithm of (3.23) has exactly the same form as (3.15). The only 
difference is that the matrix T in (3.23) is the original Toeplitz approximation matrix 
T. plus some diagonal matrix oJ. The new algorithm of (3.23) was observed to 
converge where the algorithm of (3.15) fails. It always reaches the vicinity of the true 
parameter values quickly, but requires a large number of iterations to converge to 
the exact solution (see Figure 3.4). In simulation, we observe that optimum values 
for o are inversely proportional to the condition number of R. One difficulty with 
this algorithm is that determining the optimum value for o is a problem as it is data 
dependent. Therefore, achieving optimal performance of this algorithm is dependent 
on finding a suitable value for o. We were unable to find a general expression for the 
optimum value of o in this thesis. 
4. Partitioned Toeplitz Approximation Algorithm 

An alternative to improving the condition of F in (3.15) is to block partition 
R. It may be partitioned into equally or unequally sized parts, however the number 
of row partitions must equal the number of column partitions. We now derive a 


partitioned version of (3.15). Beginning with 
Ra=r, (3.24) 


17 





Diag Pert, P=14 
e o 
a = 
“a i} 
> > 
s S 
: F 
a a 
Iterations, k Iterations, k 





Figure 3.4: Parameter tracks and general rate of convergence for the FIR 
Toeplitz approximation and Toeplitz approximation with diagonal pertu- 


bation algorithms. 


partitioning R into four parts gives 


Rn : Ry] 1 T2 


Ti + Sy : Riz ay ry 


Ra : Too + Soo} L32 T2 


Multiplying gives the following coupled equations: 


(Tir + Si1)ar + Ri2ae QT) 


(To2 + So2)a2 + Roar 


2 
and rearranging this yields 


Pr — Riza2 — Siar 


T1181 


Ty28€2 = 2 — Raa — S222. 
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(3.25) 


(3.26) 


(3.27) 


(3.28) 


But Su = Ri = Ti and S22 = Roe = To2, thus we have 


Ty ty — Ty’ (Rivae + (Ru — T1)ar) 
Tyy't2 — Tyy'(Roiar + (Roz — T22)a2) (3.29) 


ay 


a2 
and as before we let aig = Ty,'r1 and agg = Toz'f2 so that now we can write 


a = aig tai —7),'(Rizae + Rian) 
420 +ar— To (Rea + R222). (3.30) 


a2 


The final form of the partitioned Toeplitz approximation iterative algorithm is 


ar -= aigt a(*) = T(Riay” + Ruat") 
alt*?) = a2 + a”) 7 Ti (Raa\"*?) - Rnat"). (3.31) 


This partitioning algorithm was also observed to converge to the true solu- 
tion where the algorithm of (3.15) fails. Additionally, this algorithm converges to the 
true solution in about the same number of iterations as that of (3.23), and this algo- 
rithm does not require the choice of an optimal parameter value such as o in (3.23). 
It is possible to derive iterative algorithms using this same methodology beyond the 
two partitions shown in (3.25), however we did not observe a significant difference in 
performance of these algorithms over the one shown in (3.31) (see Figure 3.5). Fi- 
nally, it is also possible to combine diagonal pertubation and partitioning in a single 


algorithm, but again we did not observe any significant improvement in performance. 


B. DATA-ADAPTIVE ALGORITHMS 
Although the fixed data algorithms above do provide a solution to the least 
squares problem Ra = r, they require that a minimum amount of data be accumulated 


before the iterative algorithm may begin and a solution determined. Consequently, 
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Figure 3.5: Parameter tracks and general rate of convergence for the FIR 
partitioned Toeplitz approximation and double-partitioned Toeplitz ap- 
proximation algorithms. 

real-time processing of signals as they are available is not possible. Data-adaptive 
algorithms compute a solution based upon the inclusion of new data with each iter- 
ation. They provide a means of providing a real-time system identification solution. 
In this section we develop a data-adaptive algorithm by extending the Toeplitz ap- 
proximation technique presented in the previous section. Also, we present the well 
known RLS algorithm for comparison. 

1. Toeplitz Approximation Algorithm 


In order to avoid confusion with previous developments of the basic problem 


for the fixed data case, we define the following: 


Xn =[z(n) z(n—1) --- 2(n—M)]? (3.32) 
do 

a=|° (3.33) 
am 
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so that now y(n) = xa for a specific value of the output of an M**-order filter. For 


the general problem development, we have the following: 





z(n) z(n — 1) tee z(n — M) ao y(n) 
z(n—1 z(n—2 + 2(n-—-M-1 a n—1 
(eI) en 2 (@— m= 1) |} a} ]9— 9) ag 
r(n—P) z(n-—P-—1) --- t(n—P—M)| lay y(n — P) 
or, in vector notation we have: 
xR 
XA 
: a=y, (3.35) 
Xn-P 
and now it is easy to see that the basic problem has the form: 
xp 
: oe 
n=1 
[Xn Xn-1l --- Xp-p]] . | a=[Xn Xna-1 -*: Xnoply. (3.36) 
: ——— ee 
Xn-P : 
R 


It is important to realize how the data-adaptive algorithm can be derived from R and 


r. Using a subscript n, we can express R and r in time varying notation as: 


n 
Rn = QO xx} 
j=n—P 
= Xn_PXA_p + Xn—PaiXa—p4y + 11° + Xn-1X4_1 tXaxP 
ee alain innnnnemeaamnamemesmnall 
Rn-1 
Ry = Ra-1+%nXq (3.37) 
and 
mr, = Ss; x;y(7) 
j=n—-P 
= Xn-py(n — P) + Xn—Psiy(n — P +1) +--+ + Xn-1y(n — 1) +xny(n) 
sr” 
Tn-1 
Tn = Pn-1 + Xny(n). (3.38) 


Now that we have AR, and r,, we can form the data-adaptive Toeplitz 


approximation algorithm as in the fixed data case: 


Rna = Tp (3.39) 
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substituting R, = 7, + Sn yields 


(T, + Sy)a = Pn, (3.40) 
and rearranging gives 
Tra = Tn — Spa. (3.41) 
But S, = R, — Tp, thus 
T,a =n — (Ra —Tr)a (3.42) 


and isolating the parameter vector a we have 
a=T 'r, ~T7'Raata. (3.43) 

Finally, the iterative algorithm has the form: 
al*+!) — al 4 Tl, — T" R,al®), (3.44) 


This algorithm converges quickly to the true parameter values (see Fig 3.6). 
Notice that both algorithms begin with time n = k = 0. The data-adaptive algorithm 
begins to estimate parameter values immediately, but the fixed data algorithm must 
wait until a minimum of M +1 data points have been accumulated before it can begin 
estimating parameter values for an M‘*~order system. Although the data-adaptive 
algorithm (3.44) outperforms the fixed data algorithm of (3.31), there is a considerable 
computational tradeoff. The matricies T71 and R,, and the vector r, must be updated 
at each iteration of the data-adaptive algorithm. The matrix T7! is updated using 
the matrix inversion lemma [Ref. 10], and R, and r, are updated from (3.37) and 
(3.38), respectively. This requires approximately 5M? additional multiplications for 
each iteration of (3.44) as compared with the fixed data algorithm. 
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Figure 3.6: Parameter tracks and general rate of convergence for the FIR 
data-adaptive and fixed data Toeplitz approximation algorithms. 


2. Recursive Least Squares (RLS) Algorithm 
In order to evaluate the effectiveness of the data-adaptive Toeplitz approx- 
imation algorithm above, it must be compared against some benchmark. The RLS 
algorithm is recognized as a particularly optimal approach to the least squares FIR 


problem. A simple presentation of the update algorithm is [Ref. 5]: 
al**+t) — al*) + aR7l xne(n) (3.45) 


where a is referred to as the forgetting factor. A performance comparison between this 
algorithm and that of (3.44) reveals the astonishing fact that they perform exactly 
the same. Apppendix A shows several examples of this under different conditions of 
noise, and over- and under-modeling. On close analysis, we note that the RLS algo- 
rithm (3.45) can be expanded and rewritten in the same form as (3.44). Although 
the algorithms look and perform the same, they are in fact different. Notice that R, 
and r, are summed over the current and past data in (3.44), whereas these terms in 


the RLS algorithm are composed of only the current data values (R, = x7x, and 
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Figure 3.7: Parameter tracks and general rate of convergence for the FIR 
RLS and data-adaptive Toeplitz approximation algorithms. 

ln = Xny(n)). Interestingly, as shown in Figure 3.7, the performance of the algo- 
rithms is still the same. Unfortunately, computing a more complicated R, and r, in 
the Toeplitz Approximation data-adaptive algorithm requires an increase of approx- 
imately 2M? — M multiplications over RLS. Although this increased computation 
does not produce any improvement in performance over FIR RLS systems modeling, 


it will give improved performance over IIR RLS, as we shall see in Chapter IV. 











IV. INFINITE IMPULSE RESPONSE (IIR) 
SYSTEMS 


In this chapter we consider the infinite impulse response (IIR) system. Here 
again we develop two general categories of solution algorithms, fixed data and data- 
adaptive. The slow rate of convergence of the Gauss-Seidel iterative algorithm applies 
to IIR systems in the same way it applied to FIR systems in Chapter III. Therefore, 
in this chapter we only present one fixed data algorithm, the Toeplitz approximation 
iterative method. The data-adaptive algorithms are the recursive least squares (RLS) 
method and Toeplitz approximation method. The same discussion regarding the use 
of techniques that are Toeplitz dependent from Chapter III is applicable here as well; 


our covariance formulation of the problem precludes us from using such techniques. 


A. FIXED DATA ALGORITHMS 
1. Toeplitz Approximation Algorithm 
Here we present a Toeplitz approximation iterative algorithm [Ref. 2] of 
the same structure as in the FIR case. One major difference is that we will start with 
a partitioned version of the algorithm similar to that of (3.31). This is quite natural 


M’ represents the filter 


since we are solving for two parameter vectors. One vector, a 
input or feedforward coefficients and the other, b’, represents the filter output or 
feedback coefficients. As in Chapter III, we drop the superscripts M and N. Unless 


otherwise stated, the input coefficient vector has order M, and the output coefficient 


vector has order VV. Thus, from Chapter II, (2.18), we begin with 


Re=r (4.1) 








and partitioning R gives 


b Ty 
|-{21 «2 
Rey : Riz 2 Tr 


We now split the matricies Ry, and Rzz into a Toeplitz and a residual matrix: 


Tyy + Sy: Rye : * 
; ee ae (4.3) 
Yr 





eo eee meer er ere reese recewe 


Rey Tre + See 


Multiplication gives the following coupled equations: 


(Tyy + Syy)b + Ryza 
(Trr2 + Sez)a+ Reayb = rz (4.4) 


Ty 


and rearranging these gives 


Tyb = ry — Ryea—Syb 
Tera = Ye— Reyb— Szza. (4.5) 


But S,, = Ry — Ty, and Siz = Rez — Tez, thus 


b Ty Ty = Ty (Ryze + (Ry — Tyy)b) 


Tz re — Tz) (Reyb + (Rez — Trz)a) (4.6) 


a 


and letting bb = T~'r, and ap = T™'r, as the initial estimates of the parameter 


values, we have 


o 
tl 


bo +b — Ta"(Ryea + Ryb) 
ay +a —Tz!(Reyb + Resa). (4.7) 


~ 
il 


The final form of the fixed data IIR Toeplitz approximation iterative algorithm is: 
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Figure 4.1: Parameter tracks and general rate of convergence for the IIR 
Toeplitz approximation algorithni. 


bt) = by +bi*) — T,, (Ryza) 4 R,,b")) 
alF+)) = a + al’) —T'(R,yb4t)) ae R,ra"*)). (4.8) 


In order to evaluate the general performance of the algorithms developed 


in this chapter, the following IIR system is used for simulation: 


y(n) = 2.7600y(n — 1) — 3.8090y(n — 2) + 2.6540y(n — 3) — 0.9240y(n — 4) 
+ z(n)—0.9z(n — 1) + 0.81z(n — 2) — 0.652(n — 3) + 0.362(n — 4). (4.9) 


Figure 4.1 shows that this algorithm gives performance similar to that of the FIR 
algorithm (3.15). It converges relatively quickly when & is well-conditioned, and fails 
to converge to the true parameter values as 2 becomes ill-conditioned. 
2. Toeplitz Approximation with Diagonal Pertubation Algorithm 
As in the FIR case a significant shortcoming of (4.8) is that it fails to 
converge as the covariance matrix R becomes extremely ill-conditioned. We again 


observe that as the number of data points used to formulate R is reduced from a 
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large number (i-e., 5 to 10 times the order of the IIR filter beinz -olved for) to near 


the nuinimum number required for a unique soltion (i.e., the order of the filter}, the 
condition number of R grows quite large. Under these circumstances the algorithm 
of (4.8) does not converge to the true parameter values. We now employ the diagonal 
pertubation technique in a manner similar to that of Chapter III which improves the 
condition of R and permits a modified version of (4.8) to converge as the number of 
data points approaches the minimum required. 

If we let o, and o, be scalar multipliers and the Toeplitz approximation 
splitting of R,, and Rzz be defined as Ry, = Tyy, + Sy, and Rez = Trro + Xzz, 


respectively, then we can derive a modified version of (4.8) as follows: 
Ro=r (4.10) 


adding a diagonal matrix to Ry, and Rzz gives, 


Ry +o] : Ryz b Ty oyb 
besos ei bndig te Baek su Bibl teid Boe ae ~f=d..d 4+]... (4.11) 
Rey : Ryze +o,I} L@ rz Ora 
and substituting for Ry, and Az, we have 
Tyo t+ Sy +oyl : Ryz b Ty oyb 
pakcuv Gua bares leuntae ake es capeeeen yo es ee (4.12) 
Rey > Treo + Sez + O21 * Pe Lora 


Now, letting T,, = Ty, + o,] and T,, = Tz, + 0,1, and multiplying gives the 


following coupled equations: 


(Tyy + Sy )b+ Ryza = ry+ayb 
(Trr + Ssz)at+ Ryyb = rz+oza (4.13) 


and rearranging these gives 
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Tyb = ry, +o,b— Ryza— Syb 
Trra = Yze+o,a— Rzyb— Szza. (4.14) 


Substituting for S,, = Ry, + oyJ — Ty, and Sz, = Rrz + oz] — Tzz we have 


Tyb = ry+oyb— Ryza— (Ry + oy] — Tyy)b 
Trzpa = Yet+oza— Rzyb — (Rez + oz] — Trz)a, (4.15) 


and now we see the terms oyb and g,a cancel out: 


Tyb = ry +oyb—-oyb — Ryza— (Ry — Tyy)b 
Trpa = Y,+ora—o0,a— Rzyb — (Rzz — Trz)a. (4.16) 


After simplifying to isolate the parameter vectors, we have 


b Toye = Tyy (Ryea + (Ryy _ yy )b) 


# Yet Fs ae Tz, (Reyb + (Rez ra Tzz)a). (4.17) 


Il 


a 


Again, letting T;'ry = bo and T;;'rz = ao gives 


I 


b bo +b — T7)(Ryza + Ryyb) 


a a +a—T'(Reyb + Rrra) (4.18) 


and the final form of the Toeplitz approximation with diagonal pertubation iterative 


algorithm is: 


b+) = bo +b) — TT (Ryca® + R,,b‘*) 
a) = ap +al*) —To1(R,yb&t) + Ra). (4.19) 


The algorithm of (4.19) has exactly the same form as (4.8). Unfortunately, adding 
a diagonal value dramatically slows down the rate of convergence; however, the new 


algorithm (4.19) was observed to converge where the algorithm (4.8) fails. This 
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Figure 4.2: Parameter tracks and general rate of convergence for the IIR 
Toeplitz approximation with diagonal pertubation algorithm. 
algorithm takes a large number (i.e., several hundreds or thousands) of iterations just 
to reach the vicinity of the true parameter values, and still requires a large number of 
further iterations to converge to the correct solution. We did observe that appropriate 
values for oy and o, are inversely proportional to the condition number of R,, and 
R,z, respectively. Unfortunately, even the optimal performance of this algorithm is 
unsatisfactory (see Figure 4.2). 
3. Partitioned Toeplitz Approximation Algorithm 

In this section we block partition R,, and R,z in an attempt to improve the 
performance of (4.8). As before, the number of row partitions must equal the number 
of column partitions. Similar to previous developments we now derive a partitioned 


version of (4.8). Beginning with 


Re =r, (4.20) 





partitioning twice gives: 


(= | <= | [== | =| by Tyy | 


Ryo: | Ryo Ryzo, | Ryo bp Ty 
etimentieetiaetientnn eli eteti aie! —-— {=| -- (4.21) 
( Rey | Rey 12 | ( Rezy} Rez12 | ay Tr) 
Rey | Reyoo Rez2 | Rrx22 a2 Yr2 
Splitting the diagonal blocks into a Toeplitz and a residual matrix we have 
( Twut Sw | Ry ( Ryzy, | Ryzyy | 
Ryyor | Tyyog + Syvae Ryey, | Ryzog 
meme se nee eee x 
( Rey | Rey 12 ( Trz11 - Sez11 | Reri2 | 
Rey, | Rey Reon | Trz22 “- S2222 
| b, Ty, 
bo Tyo 
-- |=] -- |. (4.22) 
ay Tr 
a2 lr? 
Multiplying gives four coupled equations: 
(Ty + Syy11)br +f Ryy,,b2 + Rye, + Ryz 4222 = Ty, 
(Tyyoo + Syyo)b2 + Ryyo, bi + Ryegyai + Ryry922 = Pyo 
(Te211 + Sr211)a1 + Rey, bi + Rey,,b2 + Rez1282 = zy 
(Ts222 + Srx22)82 + Reyo,b1 + Reygqb2 + Reo 81 = rz (4.23) 


and rearranging yields 
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Tyyyyb1 = Ty, — Ryyygb2 — Ryey, 1 — Rye y282 — Syy),br 
Tyyqgb2 = Tyg — Ryyoybi ~ Ryzgy@1 — Ryeo9@2 — Syyogb2 ; 


Tri — Rey, 61 _ Rey,2b2 ~ Rex282 ~ Sre1181 


done a1 


Now, isolating the parameter vectors and substituting for Syy., = Ryy,,; — Tyy.1) etc., 


b, = Toe ait v1 ~ Tey yy (Rovygb2 — Ryeqy81 — Ryeyg92 — (Rey — Tyyis)br) 
bz = T yy aakv2 ~ Tey 99(Pevyo1 1 — Rye 1 — Rye,82 — (Ryyo2 — Tyy2)b2) 
a, = fie a1Tr1 ~ f Sea 11(Rey,,b1 - Rey,.b2 — Reri282 — (Reei — Tez11)a1) 
a= Te gh -T; 92( eyo, b1 — Reyogb2 — Reena: — (Re222 — Tero2)a2) (4.25) 


and letting byyg = T; wiv etc., as the initial estimate of the parameter vectors, we 


have 


by = big + by — — Ty (Reniabs — Ryz,,81 — Ryzy282 — Ryy,,b1) 
bz = bag + be — Tey an(Ryvorbs od Pyro A1 a Ryzy982 a Ryyo9b2) 
a1 = a19 + a, — eae fo ag ttey yy D2 — Rey,,b2 — Rez1282 ao Rez1181) 


a2 = A29 + ag cs Te 22(Reyo,b1 aaa Rayo, b2 iia Reo aa Rez9282). (4.26) 


The final form of the partitioned fixed data Toeplitz approximation iterative algorithm 
is: 
k k) : . ‘ 
bi) = big + bi) — Ty 11 (Ban2b2” nee )— al ) Ryyy1b)”) 
k k - ‘ 
b§ +1) bap + by? — Ty, el Banh a Rye, ~  — Ryeop )_ Ryyaob) ) 


ayt)) = aig t+ al") = Iss (Rey bf) — Rey apy Pa Rez128)” > so118\) 


= k 
af*+)) = a9 4 al") = TT ( Rib +1) -_ Reyo,b{**?) a2 Resa _ Reza"). 
(4.27 


This partitioned Toeplitz approximation algorithm was observed to con- 


verge to the true solution where the algorithm (4.8) fails, however it has the same 
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Figure 4.3: Parameter tracks and general rate of convergence for the IIR 
partitioned Toeplitz approximation algorithm. 


problem as the diagonal pertubation technique. It takes an excessive number of itera- 
tions for the algorithm to approach the vicinity of the true parameter values and even 
longer to converge exactly (see Figure 4.3). Since this is not acceptable performance 
(even marginally), we therefore turn to data-adaptive techniques in order to achieve 


the desired result. 


B. DATA-ADAPTIVE ALGORITHMS 
1. Toeplitz Approximation Algorithm 
Extending the development from equations (4.1)-(4.6) to the time-varying 


case gives the desired data-adaptive IIR Toeplitz approximation algorithm: 


bit) = b) + Ty, (Fyn a Rye, al”) ~ Ryy,b“) 
alt = alt) + Teen (Ten = Rey,b™) — Reza”) (4.28) 


where the subscript n denotes time variation. The neat and clean appearance of 


this algorithm is deceptive. Each of the terms with a subscript n must be updated 
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during each iteration. Consequently, this algorithm is somewhat computationally 
intensive. The entire algorithm is coded in matlab and is listed in Appendix D. It 
uses the matrix inversion lemma to compute the T,,,, and T;,, updates. Updating 
the rest of the time varying terms is similar to the process shown in (3.37) and 
(3.38). We observed that the algorithm gave better results for over-modeling than 
exact modeling. An example of over-modeling is where the true system output was 
generated from a 4**-order system (M = N = 4), but we are using a 6**-order model 
(M = N = 6) in the algorithm. The results for this case are shown in Figure 4.4. 
Here we see that the parameter values converge in about 15 to 20 iterations. To 
see whether this 6**-order model faithfully represents the true 4**-order system, we 
must look at the frequency response of the model versus the true system. Figure 
4.5 shows that by about 20 iterations the 6**-order model matches the frequency 
response of the true system. This performance is quite reasonable, and Appendix C 
contains parameter track and frequency response plots showing its performance under 
a variety of conditions. Another important point about this algorithm is that it does 
not use the gradient terms that appear in (2.21). In fact, when the gradient terms 
were incorporated in the algorithm, the performance was observed to be worse. 
2. Recursive Least Squares (RLS) Algorithm 

Unfortunately, there is not a well-defined algorithm that will serve as a 

benchmark against which to test the data-adaptive Toeplitz approximation algorithm. 


In this section we use, without derivation, the IIR RLS algorithm defined by[Ref. 5]: 
o (+1) — 9) 4+ aR) dae(n) (4.29) 


where @ = [b? a?]? is the parameter vector, %, = [G7 a7]? is the gradient term 
vector, and a is the forgetting factor. There are some real difficulties in keeping this 


algorithm stable [Ref. 5], and it performs best when used to exactly model a system’s 
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Data-Adaptive, M=6 N=6 
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Figure 4.4: Parameter tracks and general rate of convergence for the IIR 
data-adaptive Toeplitz approximation algorithm. 
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Figure 4.5: Frequency response of the Toeplitz approximation algorithm 
compared to the true system frequency response. 
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Figure 4.6: Parameter tracks and general rate of convergence for the IIR 
RLS algorithm. 
parameters. The performance of this algorithm, when it is stable, is shown in Figure 
4.6. However, for general system modeling where the exact order of the system is 
unknown, IIR RLS is unsatisfactory. See Figures C.10 thru C.14 of Appendix C. 
Comparing the performance of IIR RLS to that of data-adaptive IIR Toeplitz 
approximation, we observe that the increased computations required in the Toeplitz 
approximation algorithm give improved performance. The stability of the Toeplitz 
approximation algorithm is clearly evident when we consider the poles and zeros of 
the model system as it converges. In Figure 4.7 the poles (denoted by x) and zeros 
(denoted by 0) move inside the unit circle as the number of iterations goes from 9 to 
20. The fact that the poles and zeros stop moving as the algorithm continues to iter- 
ate clearly demonstrates that it is stable. Also observe that the two additional poles 
and zeros overlap at iteration k = 50, which effectively cancels their contribution to 


the frequency response. 


Similar performance is observed in the presence of noise, and the Toeplitz 





k=9 to 12 
2 k=12 to 14 


Imaginary Axis 
Imaginary Axis 





Real Axis 


k=15 to 20 


Imaginary Axis 
Imaginary Axis 





Real Axis 


Figure 4.7: Pole-zero plots for the data-adaptive IIR Toeplitz approxima- 
tion algorithm with M = N =6. 
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Freq Response at k=100 
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Figure 4.8: Parameter track and frequency response for the data-adaptive 
IIR Toeplitz approximation algorithm with M = N = 20 and 10 dB additive 
noise. 


approximation algorithm converges faster to the true frequency response when higher 
order model systems are computed. Notice how the true frequency resporse emerges 
as the algorithm begins to converge for a M = N = 20 model system in the presence 
of 10 dB additive noise (see Figure 4.9). 
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V. CONCLUSIONS 


In this thesis we applied the four concepts identified in the Introduction to 
both FIR and IIR system parameter estimation. Of these, diagonal pertubation and 
multiple partitioning are used to extend the basic fixed data algorithms of [Ref. 1] 
for FIR systems and of [Ref. 2] for IIR systems. We then used all four concepts to 
create new data-adaptive algorithms which extend the fixed data algorithms into the 
data-adaptive domain. In the following sections we draw general conclusions about 


the performance of these algorithms. 


A. FIR ALGORITHMS 

Improving the condition of the covariance matrix with increased diagonal dom- 
inance produced improved performance over the original fixed data FIR Toeplitz 
approximation algorithm. The algorithm converges fairly rapidly while using the 
minimum amount of data. Figures A.1 thru A.5 of Appendix A show the general per- 
formance of this algorithm in noise, as well as for over-modeling and under-modeling. 
The major drawback of this algorithm is the problem of determining a suitable value 
of the diagonal pertubation parameter o that will give optimal performance. 

Partitioning the covariance matrix also produced improved performance of the 
original fixed data FIR Toeplitz approximation algorithm. Using the minimum amount 
of data, the partitioned Toeplitz approximation algorithm converges almost as rapidly 
as the optimal case of the diagonal pertubation Toeplitz approximation algorithm. 
Figures A.6 thru A.9 of Appendix A show the general performance of this algorithm 
in noise, as well as for over-modeling and under-modeling. Significant features of this 


algorithm are that it does not require choosing an optimizing parameter, and the 
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partitioning reduces the computational intensity of the algorithm since the matricies 
involved are smaller. 

Extending the original fixed data FIR Toeplitz approximation algorithm to the 
data-adaptive case gives an algorithm that performs quite well, equal to RLS, however 
it is more computationally intensive than RLS. Establishing this data-adaptive FIR 
algorithm is important because it demonstrates the potential for a data-adaptive IIR 
version to work as well. Figures A.10 thru A.17 of Appendix A show the general 
performance of this algorithm and RLS under different conditions of noise, as well as 


for over-modeling and under-modeling. 


B. IIR ALGORITHMS 

Applying the same diagonal pertubation and partitioning techniques used on 
the fixed data FIR Toeplitz approximation algorithm to the IIR algorithm did not 
produce the same improvement. Although these techniques did allow the algorithm 
to converge with the desired minimum amount of data, the rate of convergence is 
unacceptable. Figures B.1 thru B.9 of Appendix B show the general performance 
of these algorithms under different conditions of data and noise, as well as for over- 
modeling and under-modeling. 

The data-adaptive IIR Toeplitz approximation algorithm succeeds in parameter 
estimation where the corresponding IIR RLS algorithms fail. It performs well under a 
wide range of noise and over-modeling conditions. We observed that even in moderate 
noise (i.e., 10 to 20 dB) over-modeling with this algorithm can be used to determine 
the true frequency response of the system being modeled. Figures C.1 thru C.14 of 
Appendix C show how the algorithms perform under different conditions of noise, 
as well as over-modeling and under-modeling. A significant consideration with this 


algorithm is the large amount of computations required for each iteration. 
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C. FUTURE WORK 

Here we have reported general observations rather than a rigorous analytical 
development on the performance of the algorithms presented in the thesis. As future 
work, the results and conclusions from this work should be put on solid analytical 
footing. Also, the algorithms presented in this thesis are developed for one-dimensonal 
real data only. The algorithms may be extended to incorporate complex and multi- 


dimensional data and related applications. 
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APPENDIX A: FIR SYSTEM SIMULATIONS 


This appendix contains plots corresponding to the FIR algorithms developed in 
Chapter III. Figure A.1 shows the performance of the unmodified Toeplitz approxi- 
mation algorithm (3.15). This plot also serves as a reference of the true system (3.8) 
frequency response. The remainder of the plcts, Figures A.2-A.17, are in the follow- 
ing general order: fixed data Toeplitz approximation with diagonal pertubation, fixed 
data partitioned Toeplitz approximation, data-adaptive Toeplitz approximation, and 


finally FIR RLS. 








6 FIR Fixed Data, P=20 
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Figure A.1: The basic fixed data FIR Toeplitz approximation algorithm. 
This plot shows the parameter tracks and frequency response of the basic Toeplitz 
approximation algorithm (3.15). The correlation matrix was formed with P = 20 
data points, and there was no additive noise in the output. The algorithm converges 
rapidly under these conditions, as can be seen by the flat parameter tracks. Also 
notice that the frequency response changes little as the number of iterations increase. 
The frequency response plot for k = 20 iterations represents the frequency response 
of the true system. 
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FIR Fixed Data, P=11 Time vs Freq Response 
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Figure A.2: Parameter tracks and frequency response of the fixed data 
Toeplitz approximation with diagonal pertubation algorithm. 

Observe that even though the minimum number of data points (P = 11) were used to 
form the correlation matrix, the algorithm converges rather quickly. The frequency 
response plot for k = 30 iterations appears identical to the response of the true system 
shown in Figure A.1. No noise was added to the output, the pertubation parameter 
was o = 1, and the true system was exactly modeled with a 10°*-order model. 
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FIR Fixed Data, P=9 


Time vs Freq Response 
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Figure A.3: Under-modeled example of the fixed data Toeplitz approxi- 
mation with diagonal pertubation algorithm. 

The 10‘*-order true system was modeled with an 8**-order model. The algorithm does 
a fairly good job preserving the frequency components of the true system, however it 
takes a much greater number of iterations to achieve this result. No noise was added 
to the output, and the pertubation parameter was o = 1. 











FIR Fixed Data, P=15 


Time vs Freq Response 
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Figure A.4: Over-modeled example of the fixed data Toeplitz approxima- 
tion with diagonal pertubation algorithm. 

The 10**-order true system was modeled with a 14**-order model. The algorithm 
converges rapidly, and by k = 30 iterations it matches the true system frequency 
response. It is important to note that this result was achieved with the minimum 
(P = 15) amount of data required. No noise was added to the output, and the 
pertubation parameter was o = 3. 





FIR Fixed Data P=21 Time vs Freq Response 
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Figure A.5: Over-modeled example with 10 dB of additive noise for the 
fixed data Toeplitz approximation with diagonal pertubation algorithm. 
The 10‘*-order true system was modeled with a 20*—order model. The additive white 
noise tends to slow down the algorithm, and distorts the frequency response it is able 
to achieve. The pertubation parameter was o = 3, and 10 dB of white noise was 
added to the output. 
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FIR Fixed Data, P=11 Time vs Freq Response 
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Figure A.6: Parameter tracks and frequency response of the fixed data 
partitioned Toeplitz approximation algorithm. 

Observe that even though the minimum number of data points (P = 11) were used to 
form the correlation matrix, the algorithm converges rather quickly. The frequency 
response plot for k = 30 iterations appears identical to the response of the true system 
shown in Figure A.1. No noise was added to the output, and the true system was 
exactly modeled with a 10**-order model. 


48 


FIR Fixed Data, P=9 Time vs Freq Response 
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Figure A.7: Under-modeled example of the fixed data partitioned Toeplitz 
approximation algorithm. 

The 10**-order true system was modeled with an 8**-order model. The algorithm does 
a fairly good job preserving the frequency components of the true system, however it 
takes a much greater number of iterations to achieve this result. No noise was added 
to the output. 
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FIR Fixed Data, P=15 Time vs Freq Response 
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Figure A.8: Over-modeled example of the fixed data partitioned Toeplitz 
approximation algorithm. 

The 10**-order true system was modeled with a 14**-order model. The algorithm 
converges rapidly, and by k = 30 iterations it matches the true system frequency 
response. It is important to note that this result was achieved with the minimum 
(P = 15) amount of data required. No noise was added to the output. 
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FIR Fixed Data, P=21 Time vs Freq Response 
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Figure A.9: Over-modeled example with 10 dB of additive noise for the 
fixed data partitioned Toeplitz approximation algorithm. 

The 10**-order true system was modeled with a 20**-order model. The additive white 
noise tends to slow down the algorithm, and distorts the frequency response it is able 
to achieve. There was 10 dB of white noise added to the output. 


51 





FIR Data-Adaptive 


Time vs Freq Response 
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Figure A.10: Parameter tracks and frequency response of the 
data-adaptive Toeplitz approximation algorithm. 

Observe that the algorithm converges to the true solution and frequency response in 
ten iterations. In this example no noise was added to the output, and the performance 
above is for an exact modeling of the 10‘*-order true system. Comparing this with 
the RLS results shown in Figure A.14 reveals that they behave exactly the same. 
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FIR Data-Adaptive Time vs Freq Res 
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Figure A.11: Under-modeled example of the data-adaptive Toeplitz ap- 
proximation algorithm. 

The 10**-order true system was modeled with an 8**-order model. The algorithm 
tries to preserve the frequency components of the true system, however it takes a 
much greater number of iterations to achieve this result. No noise was added to the 
output. 





FIR Data-Adaptive, M=14 
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Figure A.12: Over-modeled example of the data-adaptive Toeplitz approx- 
imation algorithm. 

The 10**-order true system was modeled with a 14**-order model. The algorithm 
converges rapidly, and for iterations « = 10 and higer it matches the true system 
frequency response. It is important to note that all of the additional parameters 
beyond the ten required were driven to zero by the algorithm. No noise was added 
to the output. 
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FIR Data-Adaptive, M=20 
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Figure A.13: Over-modeled example with 10 dB of additive noise for the 
data-adaptive Toeplitz approximation algorithm. 

The 10**-order true system was modeled with a 20**-order model The additive 
white noise tends to slow down the algorithm, but only slightly distorts the izequency 
response it is able to achieve. There was 10 dB of white noise added to the output. 
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Figure A.14: Parameter tracks and frequency response of the FIR RLS 
algorithm. 

Observe that the algorithm converges to the true solution and frequency response in 
ten iterations. In this example no noise was added to the output, and the performance 
above is for an exact modeling of the 10**-order true system. Comparing this with the 
data-adaptive Toeplitz approximation algorithm results shown in Figure A.10 reveals 
that they behave exactly the same. 
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Figure A.15: Under-modeled example of the FIR RLS algorithm. 
The 10‘*-order true system was modeled with an 8'*-order model. The algorithm 
tries to preserve the frequency components of the true system, however it takes a 
much greater number of iterations to achieve this result. No noise was added to the 
output. 
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FIR Data-Adaptive. M=14 


Time vs Freq Response 
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Figure A.16: Over-modeled example of the FIR RLS algorithm. 
The 10**-order true system was modeled with a 14**-order model. The algorithm 
converges rapidly, and by iteration k = 10 and higher it matches the true system 
frequency response. It is important to note that all of the additional parameters 
beyond the ten required were driven to zero by the algorithm. No noise was added 
to the output. 
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FIR Data-Adaptive. M=20 Time vs Freq Response 
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Figure A.17: Over-modeled example with 10 dB of additive noise for the 
data-adaptive Toeplitz approximation algorithm. 

The 10*-order true system was modeled with a 20**-order model. The additive 
white noise tends to slow down the algorithm, but only slightly distorts the frequency 
response it is able to achieve. There was 10 dB of white noise added to the output. 
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APPENDIX B: FIXED DATA IIR SYSTEM 
SIMULATIONS 


This appendix contains plots corresponding to the fixed data IIR algorithms de- 
veloped in Chapter IV. Figure B.1 shows the performance of the unmodified Toeplitz 
approximation algorithm (4.8). This plot also serves as a reference of the true sys- 
tem (4.9) frequency response. The remainder of the plots, Figures B.2-B.9 illustrate 
the convergence performance of the fixed data Toeplitz approximation with diagonal 
pertubation algorithm followed by the fixed data partitioned Toeplitz approximation 


algorithm. 
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IIR Fixed Data. P=100 Time vs Freq Response 
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Figure B.1: The basic fixed data IIR Toeplitz approximation algorithm. 
This plot shows the parameter tracks and frequency response of the basic Toeplitz 
approximation algorithm (4.8). The correlation matrix was formed with P = 100 
data points, and there was no additive noise in the output. The algorithm converges 
rapidly under these conditions, as can be seen by the flat parameter tracks. Also 
notice that the frequency response changes little as the number of iterations increase. 
The frequency response plot for k = 30 iterations represents the frequency response 
of the true system. 








IIR Fixed Data, P=9 
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Figure B.2: Parameter tracks and frequency response of the fixed data 
Toeplitz approximation with diagonal pertubation algorithm. 

Observe that when the minimum number of data points (P = 9) were used to form 
the correlation matrix, the algorithm converges very slowly. The frequency response 
plot for k = 50 iterations is not very close to the response of the true system shown 
in Figure B.1. No noise was added to the output, the pertubation parameters were 
oO, = 1 and o, = 10, and the true system was exactly modeled with a 4**_order model. 
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Time vs Freq Response 


IIR Fixed Data, P=7 
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Figure B.3: Under-modeled example of the fixed data Toeplitz approxi- 
mation with diagonal pertubation algorithm. 

The 4**-order true system was modeled with an 3°¢-order model. The algorithm does 
a poor job preserving the frequency components of the true system, and it takes a 
much greater number of iterations to achieve this result. No noise was added to the 
output, and the pertubation parameters were a, = 4 and a, = 50. 
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Figure B.4: Over-modeled example of the fixed data Toeplitz approxima- 
tion with diagonal] pertubation algorithm. 

The 4**-order true system was modeled with a 10*-order model. The algorithm 
converges more rapidly than the previous two examples, and by k = 100 iterations it 
nearly matches the true system frequency response. It is important to note that this 
result was achieved with the minimum (P = 21) amount of data required. No noise 
was added to the output, and the pertubation parameters were o, = 3 and o, = 100. 
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Figure B.5: Over-modeled example with 10 dB of additive noise for the 
fixed data Toeplitz approximation with diagonal pertubation algorithm. 
The 4**-order true system was modeled with a 20'*-order model. The additive white 
noise tends to slow down the algorithm, but only slightly distorts the frequency 
response it is able to achieve. The pertubation parameters were o, = 5 and a, = 100, 
and 10 dB of white noise was added to the output. 
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IIR Fixed Data, P=9 


Time vs Freq Response 
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Figure B.6: Parameter tracks and frequency response of the fixed data 
partitioned Toeplitz approximation algorithm. 

Observe that when the minimum number of data points (P = 9) were used to form 
the correlation matrix, the algorithm converges very quickly. The frequency response 
plot for k = 50 iterations is not very close to the response of the true system shown 
in Figure B.1. No noise was added to the output, and the true system was exactly 
modeled with a 4**-order model. 
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IIR Fixed Data. P=7 
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Figure B.7: Under-modeled example of the fixed data partitioned Toeplitz 
approximation algorithm. 

The 4**-order true system was modeled with an 3°¢-order model. The algorithm does 
a poor job preserving the frequency components of the true system, and it takes a 
much greater number of iterations to achieve this result. No noise was added to the 


output. 
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IIR Fixed Data, P=21 Time vs Freq Response 
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Figure B.8: Over-modeled example of the fixed data partitioned Toeplitz 
approximation algorithm. 

The 4t*-order true system was modeled with a 10*-order model. The algorithm 
converges rapidly, and by k = 20 iterations it matches the true system frequency 
response. It is important to note that this result was achieved with the minimum 
(P = 21) amount of data required. No noise was added to the output. 
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IIR Fixed Data, P=41 Time vs Freq Response 
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Figure B.9: Over-modeled example with 10 dB of additive noise for the 
fixed data partitioned Toeplitz approximation algorithm. 

The 4**-order true system was modeled with a 20‘*-order model. The additive white 
noise tends to slow down the algorithm, but only slightly distorts the frequency 
response it is able to achieve. There was 10 dB of white noise added to the output. 
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APPENDIX C: DATA-ADAPTIVE IIR 
SYSTEM SIMULATIONS 


This appendix contains plots corresponding to the data-adaptive IIR algorithms 
developed in Chapter IV. Figure C.1 shows the performance of the unmodified fired 
data Toeplitz approximation algorithm (4.8). This plot also serves as a reference of 
the true system (4.9) frequency response. The remainder of the plots, Figures C.2- 
C.4, illustrate the performance of the data-adaptive Toeplitz approximation algorithm 


followed by the IIR RLS algorithm. 
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IIR Fixed Data. P=100 Time vs Freq Response 
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Figure C.1: The basic fixed data IIR Toeplitz approximation algorithm. 
This plot shows the parameter tracks and frequency response of the basic Toeplitz 
approximation algorithm (4.8). The correlation matrix was formed with P = 100 
data points, and there was no additive noise in the output. The algorithm converges 
rapidly under these conditions, as can be seen by the flat parameter tracks. Also 
notice tnat the frequency response changes little as the number of iterations increase. 
The frequency response plot for k = 30 iterations represents the frequency response 
of the true system. 








Data-Adaptive, M=4 N=4 
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Figure C.2: Parameter tracks and frequency response of the data-adaptive 
Toeplitz approximation algorithm. 

Observe that the algorithm converges to the true solution and frequency response 
in about 60 iterations. In this example no noise was added to the output, and the 
performance above is for an exact modeling of the 4'*-order true system. Comparing 
this with the RLS results shown in Figure B.18 reveals that they do not behave the 
same as in the FIR case. 
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Figure C.3: Pole-zero plots for the data-adaptive Toeplitz approximation 
algorithm. 

Observe that as the algorithm converges to the true solution, the poles (x’s) move 
inside the unit circle. After about 50 iterations they stop moving, and their location 
corresponds to the frequency response shown in Figure B.10. Finally, we see that 
there is almost no movement as the algorithm iterates from k = 70 to 100 iterations. 
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Figure C.4: Under-modeled example of the data-adaptive Toeplitz approx- 


imation algorithm. 
The 4'*-order true system was modeled with a 3"¢-order model. The algorithm tries 


to preserve the frequency components of the true system, but with only three poles 
and zeros it is unable to preserve both of the main components. No noise was added 


to the output. 
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Figure C.5: Pole-zero plots for the under-modeled data-adaptive Toeplitz 
approximation algorithm. 

Observe that as the algorithm converges, the poles (x’s) still move inside the unit 
circle. After about 50 iterations their location corresponds to the frequency response 
shown in Figure B.12. Finally, we see that the pole-zero pair on the real axis tend to 
cancel each other, which explains why there is only one primary frequency component 
in Figure B.12. 


75 








Data-Adaptive, M=10 N=10 
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Figure C.6: Over-modeled example of the data-adaptive Toeplitz approx- 


imation algorithm. 


The 4t*-order true system was modeled with a 10‘*-order model. The algorithm 
converges rapidly, and by k = 50 iterations it matches the true system frequency 
response. It is important to note that all of the additional pole zero pairs beyond the 
four required were driven to cancel each other by the algorithm. No noise was added 


to the output. 
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Figure C.7: Pole-zero plots for the over-modeled data-adaptive Toeplitz 
approximation algorithm. 

Observe that as the algorithm converges, the poles (x’s) are all inside the unit circle. 
After about 50 iterations their location corresponds to the frequency response shown 
ir. Figure B.14. Finally, we see that the pole-zero pairs different from those shown in 
Figure B.11 are overlapping and tend to cancel each other, which explains why both 
of the dominant frequency components in Figure B.10 are present. 
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Data-Adaptive, M=20 N=20 Time vs Freq Response 
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Figure C.8: Over-modeled example with 10 dB of additive noise for the 
data-adaptive Toeplitz approximation algorithm. 

The 4**-order true system was modeled with a 20**-order model. The additive white 
noise teids to slow down the algorithm, but only slightly distorts the frequency 
response it is able to achieve. There was 10 dB of white noise added to the output. 
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Figure C.9: Pole-zero plots for the over-modeled data-adaptive Toeplitz 
approximation algorithm with 10dB additive noise. 

Observe that as the algorithm converges, the poles (x’s) are all inside the unit circle. 
After about 100 iterations their location corresponds to the frequency response shown 
in Figure B.16. Again, we see that the pole-zero pairs different from those shown in 
Figure B.11 are almost overlapping and tend to cancel each other, which explains 
why both of the dominant frequency components in Figure B.10 are present. 
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Data-Adaptive, M=4 N=4 Time vs Freq Response 
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Figure C.10: Parameter tracks and frequency response of the IIR RLS 
algorithm. 

Observe that the algorithm converges to the true solution and frequency response in 
ten iterations. In this example no noise was added to the output, and the performance 
above is for an exact modeling of the 4**-order true system. Comparing this with the 
data-adaptive Toeplitz approximation algorithm results shown in Figure B.10 reveals 
that they behave differently. 
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Figure C.11: Pole-zero plots for the IIR RLS algorithm. 
Observe that as the algorithm converges to the true solution, the poles (x’s) move 
inside the unit circle. After about 15 iterations they stop moving, and their location 
corresponds to the frequency response shown in Figure B.18. Finally, we see that 
there is almost no movement as the algorithm iterates from k = 18 to 50 iterations, 
which indicates that the algorithm has stabilized. 
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Figure C.12: Under-modeled example of the IIR RLS algorithm. 
The 4'*-order true system was modeled with an 3°?-order model. The algorithm 
tries to preserve the frequency components of the true system, however it becomes 
unstable after about 50 iterations. No noise was added to the output. 
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Figure C.13: Pole-zero plots for the under-modeled IIR RLS algorithm. 
Observe that as the algorithm converges, the poles (x’s) are still outside the unit 
circle. After about 50 iterations their location reflects the fact that the algorithm is 
unstable, and it fails to converge. 
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Figure C.14: Over-modeled example of the IIR RLS algorithm. 
The 4'*-order true system was modeled with a 10**-order model. The algorithm also 
fails to converge, as can be seen after about 29 iterations in the parameter track plot. 
It is worth noting that for a brief number of iterations, from about k = 10 to 27, the 
algorithm does produce the true system frequency response. Unfortunately, this is 
not a stable result of the algorithm. No noise was added to the output. 
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APPENDIX D: MATLAB CODE FOR IIR 
ADAPTIVE ALGORITHM 


% This is the IIR ADAPTIVE algorithm which uses the Toeplitx Approximation 


% splitting R=T+S. 


N = input(’How many poles in the system model would you like, N =? ’); 
M = input(’How many zeros in the system model would you like, M=7? ’); 
k = input(’Enter the number of iterations (i.e. k = 100). ’); 

noise = input(’How much noise would you like (i.e., noise = .1) ? ’); 


maxMN = max(M+1,N); 


MN = M +N; 
s = .0000001; % A small amount used to start the algorithm below. 
L = 1; 


P=1; % Sart with one data point. 


% Now generate the "true" system data. 


rand(’normal’); 


rand(’ seed’ ,0); 


a = [1,-2.76,3.809,-2.654,.924]; % These are the true system 
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b = [1,-.9,.81,-.65, .36]; % coefficients. 

x = [zeros(M+1,1);rand(P + k,1)]; % The input data, white noise. 

y = filter(b,a,x); % The output data. a 
x = flipud(x); : % Reorder the input for convenience. 
rand(’seed’ ,5) ; % A separate noise generator. 

y = y + sqrt(noise) * [zeros(M+1,1);rand(P+k,1)]; % Add noise to output. 

y = flipud(y); 

y = [y(i:P+k,1);zeros(maxMN+1,1)]; % The additional zeros are used in the 

x = ([x(1:P+k,1);zeros(maxMN+1,1)]; % computation of Z below, read it. 


% Now I will construct data matrix, Z 


Z = zeros(P,MN + 1); 


for i = 1:P 

Z(i,:) = [y(itk+2:N+k+1+i)? ,x(itk+1:i+k+M+1) J; 

% The counting index 

% ensures that the algorithm starts with zero data points. 


end 


% Now generate the data correlation matriz, R, 


% and partition it into 4 parts given by M and N. 


RZ’ * Z; 


Ryy = R(1:N,1:N); 








Ryx = R(1:N,N+1:M+N+1) ; 
Rxy = Ryx’; 


R(N+1:M+N+1,N+1:M+N+1); 


i 


% Now generate the Toeplitz approximation to Ryy and Rxx. 


tyy = zeros(N,1); 
txx = zeros(M+1,1); 
for i= 1:N 


~ 


tyy(i) = mean(diag(Ryy,i-1)); 


end 
for i = 1:M+1 
txx(i) = mean(diag(Rxx,i-1)); 


end 


Tyy = toeplitz(tyy) + s * eye(N); % s is added so that 


Txx = toeplitz(txx) + s * eye(M+1); % these may be inverted. 


Tiy = inv(Tyy); 


Tix = inv(Txx); 


% Now generate the cross-correlation vector, r, 


% and partition it into two parts given by M and N. 
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r= Z?’ * y(2+k:P+k+1) ; 


ry = 


r(1:N); 


4 
* 
" 


r(N+1:M+N+1) ; 


% Now compute the initial estimates of the parameter vectors. 


oC 
i 


Tiy * ry; 


a= Tig *-2s; 

Qk = zeros(k,M+N+1); % Holds the results of each iteration. 
for i= Ack; % This is the algorithm!!! 

Xn = y(k+3-i : k+N+2-i); % get new output data 

Xm = x(k+2-i : k+2+M-i); % get new input data 


Yn = y(k+2-i); % get next output data point 


% Use the matrix inverse lemma to compute the next 


% values for Tiy and Tix, note that L=1. 


Tiy = (Tiy - ( Tiy * Xn * Xn’ * Tiy) / (L+Xn’*Tiy*Xn))/L; 
Tix = (Tix - ( Tix * Xm * Xm’ * Tix) / (L+Xm’*Tix*Xm))/L; 


% Update all of the rest of the time-varying quantities. 


Ryy = Ryy + Xn * Xn’; 
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Ryx 
Rxy 


Rxx 


% Co 


o 
ul 


p» 
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Qk(i 


end; 


+ Xn * Xm’; 


W 
a 
Fa] 


= Rxy + Xm * Xn’; 
= Rxx + Xm * Xm’; 
ry + Xn * Yn; 


rx + Xm * Yn; 


mpute the next value of the parameter vectors. 


b + Tiy * (ry - Ryx * a 





Ryy * b); 


a+ Tix * (rx - Rxy * b - Rxx * a); 


,:) = [b’,a’]; 
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